Phase-field crystal study of grain-boundary premelting 
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We study the phenomenon of grain-boundary premelting for temperatures below the melting 
point in the phase-field crystal model of a pure material with hexagonal ordering in two dimensions. 
We investigate the structures of symmetric tilt boundaries as a function of misorientation 6 for 
two different inclinations and compute in the grand canonical ensemble the "disjoining potential" 
V(w) that describes the fundamental interaction between crystal-melt interfaces as a function of the 
premelted layer width w, which is defined here in terms of the excess mass of the grain boundary 
via a Gibbs construction. The results reveal qualitatively different behaviors for high-angle grain 
boundaries that are uniformly wetted, with w diverging logarithmically as the melting point is 
approached from below, and low-angle boundaries that are punctuated by liquid pools surrounding 
dislocations, separated by solid bridges. The latter persist over a superheated range of temperature. 
This qualitative difference between high- and low- angle boundaries is reflected in the w-dependence 
of the disjoining potential that is purely repulsive (V'(w) < for all w) for misorientations larger 
than a critical angle 8 C , but switches from repulsive at small w to attractive at large w for 9 < 9 C . 
In the latter case, V(w) has a minimum that corresponds to a premelted boundary of finite width 
at the melting point. Furthermore, we find that the standard wetting condition 7 g b(# c ) = 27 s ; gives 
a much too low estimate of 8 C when a low-temperature value of the grain boundary energy y s b is 
used. In contrast, a reasonable lower-bound estimate can be obtained if 7 g i, is extrapolated to the 
melting point, taking into account both the elastic softening of the material at high homologous 
temperature and local melting around dislocations. 

PACS numbers: 61.72.Mm, 64.70.D-, 81.16.Rf, 81.30.Fb 
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I. INTRODUCTION AND SUMMARY 

The presence of liquid films at grain boundaries for 
temperatures below the melting point can alter macro- 
scopic properties of polycrystalline materials and dra- 
matically reduce resistance to shear stresses. The latter 
can lead to catastrophic material failure as exemplified 
by hot cracking during high-temperature processing of 
metallic alloys^ 2 -. While there is indirect experimental 
evidence for the occurrence of grain-boundary premelt- 
ing in both pure materials 3 -^ and alloys^, it is inher- 
ently difficult to image and to measure thermodynamic 
properties of nanometer-width liquid films. One excep- 
tion is optical microscopy of colloidal crystals, which has 
produced striking "atomistic" -scale images of premelted 
grain boundaries^. Even in this case, however, the lack 
of precise control of grain geometry and external condi- 
tions makes it hard to determine the fundamental nature 
of the premelting transition. 

Molecular dynamics (MD) simulations provide in prin- 
ciple a powerful alternative to experiment for studying 
grain-boundary premelting. MD studies using Lennard- 
Jones^^ and interatomic potentials for metal s 10 ' 11 ^ 2 
and semiconductors such as silicon— have reported ev- 
idence for disordered layers at grain boundaries at dif- 
ferent temperatures below^^ i 10 ! 11 : 13 and abovei^ the 
melting point. In addition, such layers have been re- 
ported to exhibit fluid-like properties in a MD study 
of grain-boundary shearing in a Lennard-Jones system 
where the shear modulus decreased sharply below the 



melting point^. The large fluctuations inherent in MD 
simulations, however, make it generally hard to compute 
precisely the thermodynamic properties of grain bound- 
aries at high homologous temperature and to quantify 
the interaction between crystal-melt interfaces. 

A remarkable study of grain-boundary premelting was 
carried out by Kikuchi and Cahnii using a lattice gas 
model and a cluster variation approximation for the eval- 
uation of its thermodynamic properties. Their results 
were later corroborated by Monte Carlo simulations of 
the same models. They indeed found a liquid-like layer 
at the grain boundary for temperatures well below the 
melting point. The width of this layer diverges logarith- 
mically when the melting point is approached. While this 
study gave valuable insights, it did not yield a complete 
picture of grain-boundary premelting since the construc- 
tion of the lattice-gas model leads to numerous geomet- 
rical constraints, such that only a single misorientation 
could be investigated. 

From a basic thermodynamic viewpoint, grain- 
boundary premelting is governed by the balance between 
bulk and interfacial free energies. While the difference 
in bulk free energies per unit volume between solid and 
liquid^, AG(T) = G S (T) - Gi(T), always favors a crys- 
talline state below the melting point, the interfacial free 
energy favors the formation of a liquid layer for wetting 
conditions. The total excess free energy (per unit area of 
grain boundary) that reflects both contributions can be 



written in the form 17 



A. Sharp- and diffuse-interface theories 



G CXC ( W ,T) = AG(T)w + 2 lsl + V(w), 



(1) 



where w denotes the liquid layer width and the last two 
terms on the right-hand-side represent the interfacial free 
energy. The latter must reduce to twice the excess free 
energy of the solid- liquid interface, 27,,;, when the two 
solid-liquid interfaces are well separated, but generally 
contains an additional contribution V(w) when their sep- 
aration becomes comparable to the intrinsic nanometer- 
width 5 of an isolated solid-liquid interface. This ad- 
ditional contribution, referred to hereafter as the "dis- 
joining potential," represents the interaction due to the 
overlap of two solid-liquid interfaces, which drives the 
formation of a liquid layer under wetting conditions, or 
conversely joins two crystals for non-wetting conditions. 
Its derivative V'(w) is directly analogous to the disjoining 
pressure used in the physics of thin liquid films 18 . 

In this paper, we study grain-boundary premclt- 
ing using the phase-field crystal (PFC) modeling 
approac h 19 ' 20 ! 21 ' 22 ! 23 ' 24 inspired from classical density- 
functional theory— In the present context, this mean- 
field approach has the advantage of resolving the atomic- 
scale density-wave structure of a polycrystalline mate- 
rial while, at the same time, averaging out fluctuations. 
Therefore, it is ideally suited for computing quantita- 
tively the disjoining potential and for elucidating its 
relationship to atomic grain-boundary structure. In a 
recent study, Berry et al. observed melting at grain 
boundaries in three-dimensional phase-field crystal sim- 
ulations for bcc ordering 26 , thereby suggesting the use- 
fulness of this method for investigating fundamental as- 
pects of this phenomenon. Thermodynamic properties 
of premelted grain boundaries, however, were not stud- 
ied in detail in this work. The present work focuses 
on the quantitative study of premelting in the gen- 
eral framework of Eq. {T]) with the appropriate choice 
of thermodynamic ensemble for the phase-field crystal 
model, for two-dimensional crystals with hexagonal sym- 
metry. We compute explicitly the dependence of the 
disjoining potential on layer width and determine wet- 
ting conditions as a function of grain-boundary orien- 
tation. This allows us to make contact with sharp-i 8 - 
and diffuse-interface 1 - ' 27 ' 28 theories of interfacial premelt- 
ing. Like MD studies with truncated short-range in- 
teratomic potentials^' 9 - ' 10 ' 11 ' 12 ' 13 , these theories neglect 
the effects of long-range dispersion forces considered in 
statistical theories of grain-boundary melting^ 9 - and in 
theoretical 3 ^ and experimental 3 - 1 - studies of intergranu- 
lar phases in ceramic materials. These forces are also 
neglected in the present phase-field crystal study that fo- 
cuses on the structural component of the disjoining po- 
tential due to partial crystal ordering within premelted 
layers. 



In the simplest theory, a "wet" grain boundary is 
modeled as a thin layer of liquid sandwiched between 
two solid-liquid boundaries, assumed to be sharp and 
straight. If only short-range forces are present, an ex- 
ponential interaction between the interfaces is expected 
for large film thickness 1 ^. This suggests to write 1 ^ 



V(w) = A7exp(--J , 



(2) 



where the prefactor A7 = 7° b — 27 s j guarantees that 
the total interfacial free energy V(w) + 2j s i in Eq. (Q]) 
reduces to the energy 7° b of a "dry" grain boundary in 
the limit of vanishing liquid layer width. Minimization of 
the total excess free energy in Eq. |T]) with respect to w, 
with V(w) defined by Eq. @, predicts that, for A7 > 
0, the liquid layer width vanishes for temperatures less 
than a "bridging temperature" T&, defined by AG(Xb) = 
— A7/5, and increases smoothly as 

w{T) = -8ln(-AG(T)S/Aj) for T b < T < T m , (3) 

ultimately diverging as the melting temperature T m is 
approached from below. For A7 < 0, in contrast, bound- 
aries remain completely dry (w(T) — 0) for all T up to a 
maximum temperature T* defined by 



AG(T* 



-Aj/5, 



(4) 



and are in metastable equilibrium with respect to the 
liquid in the superheated range T m < T < T* . 

The grain-boundary energy is generally defined as the 
total excess free energy of the boundary with respect to 
the solid, or 7 gfc (T) = G exc {w(T), T) here. Therefore, 
in the wetting case (A7 > 0), this energy is constant 
and simply equal to 7° 6 for T < T^, consistent with the 
requirement that 7 g fc(0) = 7° b , but decreases for T > Tb 
until reaching 2j s i at the melting point where w diverges. 
Substituting Eq. © into Eq. Q gives 

7gb (T) - 2 lsl = -SAG(T) [1 + In (-AG(T)5/A 7 )] , (5) 

for Tb < T < T m , which has limits A7 and zero at T = Tb 
and T = T m , respectively. In contrast, for non- wetting 
conditions (A7 < 0), this theory predicts that the grain- 
boundary energy retains its dry value for all tempera- 
tures: 7g6 (T) = 7°, for T < T*. 

Phase-field theories of interfacial premelting 27,28 where 
solid-liquid interfaces are inherently spatially diffuse have 
yielded predictions that are in part consistent with the 
above picture, but also point to the possibility of more 
complex premelting behaviors. The most detailed stud- 
ies have been carried out in models where the crystal 
orientation is represented by a scalar field coupled to the 
standard scalar phase field that measures the local crystal 
disorder. For wetting conditions, those models predict ei- 
ther a smooth increase of w with temperature below T m , 
qualitatively similar to the behavior predicted by Eq. ^ , 



or the existence of first-order transitions between grain- 
boundary states of different width s 27 ' 28 , in analogy with 
the theory of critical-point wetting^ 2 -. These predictions, 
however, depend generally on the choice of phcnomeno- 
logical thermodynamic functions and parameters in those 
models that cannot be derived directly from microscopic 
physics, and thus are hard to relate to real systems. The 
phase-field crystal model, in contrast, has the advantage 
of removing much of the arbitrariness inherent in conven- 
tional phase-field theories, ft explicitly describes the dis- 
location structure of grain boundaries and is formulated 
in terms of physical quantities, such as the liquid struc- 
ture factor, that can be either measured experimentally 
or computed using MD simulations. Hence, this model 
can in principle make quantitative predictions that can 
be compared to both experiments or MD simulations as 
demonstrated recently for isolated solid-liquid interfaces 
in a bcc syste m 33 ' 34 . 



B. Disjoining potential and layer width definitions 
in the phase-field crystal model 

Before summarizing our main results, some thermody- 
namic considerations relevant for the present phase-field 
crystal study are worthy of brief mention. First, while 
premelting for pure materials has traditionally been dis- 
cussed in the Gibbs ensemble of Eq. ([TJ) with constant 
T, pressure p, and particle number N, the choice of the 
grand canonical ensemble with constant T, chemical po- 
tential fi, and volume V, is more suited for the PFC 
model where the Helmholtz free energy is a function of 
the density and simulations are carried out at fixed V. 
Thus, the disjoining potential is defined here in terms 
of the excess of the grand potential in complete analogy 
with Eq. fT}, i.e. it represents the total interfacial con- 
tribution of this excess minus its asymptotic value for 
well-separated interfaces equal to 2j s i. For reasons de- 
tailed below, it is simpler to study premelting as a func- 
tion of /i rather than T in the PFC model. Both are 
intensive variables, and the results are expected to be 
equivalent, fn particular, the departure of the chemical 
potential from its equilibrium value /z oq — /j, in the grand 
canonical ensemble is analogous to the departure of the 
temperature from the melting point T — T m in the Gibbs 
ensemble, with the solid (liquid) being stable for negative 
(positive) values of both quantities in their respective en- 
sembles. For convenience, even though we work in the 
grand canonical ensemble with \i as control parameter, 
we often refer interchangeably hereafter to temperature 
and chemical potential to facilitate the comparison of our 
results to previous theories and experiments. 

Second, we define the liquid layer width using a Gibbs 
construction. We first determine the excess mass car- 
ried by the grain boundary, which is simply the total 
mass of the bicrystal system with a grain boundary at 
fixed jU minus the mass of a single crystal occupying the 
same volume at the same (J,. The film thickness w is then 



defined by equating this excess mass to the product of 
w and the difference between solid and liquid densities. 
The advantage of this thermodynamic approach is that 
it gives a precise definition of w that remains applicable 
even when the liquid layer is not spatially uniform along 
the grain boundary, as is the case here for small misorien- 
tations. This definition of course reduces to the standard 
definition of the layer width in the limit where the liquid 
layer width is much larger than the intrinsic solid-liquid 
interface width (w > 5). 



C. Main results 

Let us now summarize our main results as they re- 
late to the theories reviewed above. The structure and 
properties of symmetric tilt boundaries were studied as 
a function of misorientation 9 for two different inclina- 
tions where the symmetry axis (from which each crystal 
is rotated by ±0/2) is parallel (<fi = 0) or at a 30° an- 
gle (4> = 30°) to any of the six equivalent close-packed 
directions of the hexagonal crystal. 

We find that high-angle boundaries behave essentially 
as predicted by the sharp-interface theory. They are dry 
well below the melting point and become uniformly wet- 
ted with a liquid layer of roughly constant width along 
the boundary. The latter diverges logarithmically as the 
melting point is approached from below, consistent with a 
disjoining potential that is reasonably well approximated 
by the simple exponential form of Eq. ([2]). fn contrast, 
the behavior of low-angle boundaries in the PFC model 
is not correctly predicted by the sharp-interface theory, 
both qualitatively and quantitatively. The main qual- 
itative difference is that grain boundaries in the PFC 
model do not remain dry with zero width as predicted 
by this theory. They exhibit local melting around dis- 
locations, as previously seen in Rcf. 26 , and the resulting 
liquid pools cause w to increase smoothly with tempera- 
ture (i.e., fi eq — fj, in the grand canonical phase-field crys- 
tal simulations) , although w remains finite at the melting 
point and into the superheated range (/x eg — /U > 0) where 
these boundaries are metastable. At a more quantita- 
tive level, dislocation-induced premelting contributes to 
the reduction in the grain-boundary energy from its low- 
temperature value, which can be larger than 2 , y s / even 
for small misorientations, to a value less than 2%i near 
the melting point. The other factor contributing to this 
reduction is the elastic softening of the material at the 
melting point discussed below. 

Dislocation-induced premelting of low-angle bound- 
aries is reflected in the ^-dependence of the disjoining 
potential V(w) that exhibits a minimum at a finite width 
w = w m , which corresponds to the equilibrium layer 
width at the melting point. Therefore, this potential is 
repulsive for w < w m and attractive for w > w m . In 
contrast, it is predicted to be attractive for all w in the 
sharp-interface theory. The high- and low-angle regimes 
can be formally distinguished by defining a critical mis- 



orientation C such that, for > C , the disjoining poten- 
tial is purely repulsive for all w, and, for < C , exhibits 
a minimum with short-range repulsion and long-range 
attraction. Our results suggest that the transition be- 
tween these two regimes is smooth, with the equilibrium 
layer width at the melting point diverging in the limit 
where 9 approaches 9 C from below, although the nature 
of this divergence is hard to pinpoint precisely. It should 
be emphasized that the transition does not correspond 
to a sharp transition in the geometry of the grain bound- 
ary. Rather, the critical angle falls into a range where 
the geometry of the grain boundary is somewhere in be- 
tween the two extremes described above. Namely, when 
the melting point is approached from below for 9 slightly 
above or below C , the grain boundary consists of liquid 
pools separated by "bridges" , but the distance between 
the dislocations is comparable to the pool diameter, so 
that the pools start to overlap and the material of the 
bridges is no longer fully solid. 

While the main results described so far are ostensibly 
independent of inclination, we find some additional in- 
dependent features that require refinement of the above 
picture. For the <j> = 30° inclination, which has the 
simplest behavior, the liquid pools were always found 
to be centered around isolated dislocations for low-angle 
boundaries, as seen qualitatively in Ref^, and to merge 
progressively to form a uniform film with increasing mis- 
orientation. In contrast, for the <j) — inclination, 
which was investigated here in greater detail, discontin- 
uous structural transitions were seen between different 
grain-boundary states. In one state, each dislocation is 
surrounded by its own liquid pool. In the other state, 
two dislocations combine to share a common liquid pool. 
The structural transition between these two states only 
occurs above a small misorientation well below 9 C . The 
transition first occurs in the overheated states above the 
melting point, and shifts to lower values of /u e q — A* with 
increasing misorientation. Furthermore, the transition is 
hysteretic, such that two grain-boundary states with dif- 
ferent liquid-pool structures can coexist over a certain 
range of /i. The jump in w at the transitions between 
different states, which measures effectively the change in 
liquid fraction associated with the pairing of liquid pools, 
is small. Hence, it cannot be ruled out that these transi- 
tions would be smeared by fluctuations and not directly 
observable as sharp transitions in a real system. 

The most relevant aspect of our results for experiment 
is the quantitative prediction of the critical wetting an- 
gle 9 C above which the liquid layer width diverges at 
the melting point. In the sharp-interface theory, this 
angle is predicted by the standard wetting condition 
A7(0) = 7° 6 (0) - 27 s/ = 0, where 7° b (0) is taken to 
be the completely dry grain-boundary energy far below 
the melting point. As such, this condition predicts a 
value of 9 C that is much smaller than observed in the 
phase-field crystal simulations. This failure is due to 
the fact that the sharp-interface theory predicts that the 
grain-boundary energy is constant for nonwetting condi- 



tions. As noted earlier, this energy is reduced by both 
dislocation-induced premelting and the elastic softening 
of the material at high homologous temperature. There- 
fore, one would expect a better estimate of 9 C to be ob- 
tained by comparing 2"/ s i to a value of the grain-boundary 
energy at the melting point, 7™ b (0), which is generally 
much lower than 7° b (0) as found in a MD study of a tilt 
boundary in pure Cv*2£. Note that 7X(0) — * 2j s i when 9 
approaches 9 C from below, and 7I? = 2^ s i for all 9 larger 
than 9 C . 

Of course, the precise determination of 7?X(0) gener- 
ally requires a complete solution of the problem since it 
depends on the structural details of the premelted grain- 
boundary structure. A somewhat better prediction of 
9 C can nonetheless be obtained by an estimation of the 
grain-boundary energy at the melting point that takes 
into account the bulk elastic softening of the material 
and melting around dislocations. As in previous PFC 
studie s 19 ' 20 , we find that for low-angle grain boundaries 
7 g fc is well described by the Read-Shockley law2£. The 
physical parameters entering this law are the shear mod- 
ulus G and the dislocation core radius r$. Hence we have, 
for small angles, 7 g {, ~ 7rs(0, G, ro). Elastic softening 
and dislocation premelting are reflected in the tempera- 
ture dependence of these quantities. The shear modulus, 
which can be calculated analytically in the PFC model, 
has large variations: denoting by Go and G m its values 
at zero temperature and at the melting point, respec- 
tively, we find typically Go/G m ~ 3. The variation of 
the core radius is obtained by fitting our simulation data 
with the Read-Shockley law and describes phenomeno- 
logically the dislocation premelting. We observe an in- 
crease of the core radius at the melting point by about 
40% with respect to its zero-temperature value. 

It should be noted that, whereas data for the vari- 
ation of the elastic constants are readily available, the 
variation of the core radius is a result of the premelt- 
ing around dislocations and hence difficult to quantify 
in experiments or MD simulations. This suggests that 
it is useful to consider two successive approximations to 
improve the estimate of the critical wetting angle. If 
only the elastic softening is included, we can exploit the 
fact that in the Read-Shockley law the grain-boundary 
energy is simply proportional to the shear modulus. 
Therefore, we have 7^(0) w 7° b (0)G ro /G o , and thus 
the modified wetting condition 7° b (0 c )G m /Go ~ 2jsi- 
If, in addition, the variation of the core radius is in- 
cluded, the estimate for the grain-boundary energy be- 
comes 7^(0) «7° & 7 R s(0,G m ,r o (T m ))/7 RS (0, Go, r (0)), 
where ro(0) and ro(T m ) are the values of the core radius 
at zero temperature and the melting point, respectively. 
Inserting the explicit expression of the Read-Shockley law 
yields the improved wetting condition 

where a is the lattice constant of the hexagonal crys- 



tal, and a — s/3/2 is the distance between close-packed 
planes in the hexagonal structure, expressed in units of 
the lattice spacing. Concretely, for the = inclina- 
tion, the phase-field crystal simulations yield 9 C w 14°. 
The standard wetting condition 7° fc = 2^f s i predicts a 
completely erroneous value of 9 C of about 2° . The condi- 
tion including only the elastic softening predicts ff c w 6°, 
whereas the condition of Eq. © including both effects 
yields 9 C « 10°. The remaining discrepancy with the 
value from simulations reflects the fact that the Read- 
Shockley law is no longer valid when liquid pools start 
to overlap and the structure of the grain boundary is no 
longer well described by an array of isolated dislocations, 
which is precisely the range where 9 ~ 9 C . While, there- 
fore, even the best estimate is still of limited accuracy, 
it sheds light on several physical effects determining the 
critical wetting angle that have not been previously ap- 
preciated. 

Finally, the failure of the sharp- interface theory to pre- 
dict the critical wetting angle obviously makes this theory 
inadequate in predicting the superheated range of tem- 
perature for low- angle boundaries. For = 0, this theory 
predicts that only boundaries with 9 less than 9 C « 2° 
are superheated, while boundaries in the PFC model can 
be superheated for 9 up to 9 C i=s 14°, with this range 
vanishing as 9 — > 9 C . If the computed values for the 
grain-boundary energy at the melting point arc used in- 
stead of the low-temperature values, the sharp-interface 
prediction yields the right order of magnitude for the su- 
perheated range for angles close to 9 Cl but largely over- 
estimates this range for low-angle grain boundaries. 

The remainder of the paper is organized as follows. In 
Sec. II, we briefly review the phase-field crystal model 
and describe our numerical methods. In Sec. Ill, we out- 
line the procedure for obtaining the liquid layer width, 
the grain-boundary energy, and the disjoining potential 
from our simulation data. Our results are presented in 
Sec. IV, and discussed in Sec. V. Finally, concluding re- 
marks and future prospects are given in Sec. VI. 



II. PHASE-FIELD CRYSTAL MODEL 

A. Basic equations and properties 

We consider the simplest PFC model defined by the 
dimensionless free energy functiona l 19 ' 20 



^ = /^{!H + (v 2 + i) 2 ]^ + ^ 4 } 



(7) 



which is a transposition to crystalline solids of the Swift- 
Hohenberg model of pattern formation 3 -!. Furthermore, 
we define the dimensionless chemical potential 



/' 
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-eip + (V 2 + 1) 2 t/> + tp 3 



(8) 



The dimensionless functional in Eq. ([7]) can be obtained 
by a suitable rescaling of a dimensional free-energy func- 



tional which, in turn, can be related to classical density- 
functional theory. Since these transformations have been 
discussed in detail elsewher e 19 ' 20 ' 24 ' 33 ' 34 , they do not 
need to be repeated here. 

The phase diagram of this model has also been dis- 
cussed previously 2 ^. However, since a precise character- 
ization of the bulk phases is important for the present 
work, we resume here the main steps that are necessary 
to obtain the properties which are needed in the subse- 
quent developments. To construct the phase diagram, we 
calculate separately the free-energy density (free energy 
per unit surface in two dimensions) as a function of the 
mean density ip in the solid, denoted by f s (ip), and in the 
liquid, /z(V>), using Eq. (|7|). Since the density is uniform 
in the liquid, fi(ip) is obtained directly from Eq. ([7]), 



/i = -(e-l)£ + £. 

•" v ' 2 4 



(9) 



It is possible to obtain an analytical expression for f s in 
the one-mode approximation, in which only the contri- 
bution of the principal reciprocal- lattice vectors is taken 
into account. Then, the density for the two-dimensional 
hexagonal structure can be written^ as 



tps(x,y) = ip + A t 



, , 1 iy\ 1 { 2 <iy 

cos(qx) cos —= cos —= 



(10) 

This solution ansatz is inserted into the free energy, 
Eq. ([Jj. Integrating over a unit cell and minimizing the 
free energy with respect to A t and q leads to 



*-!N 



15e - 36ip : 



(11) 



where the ± signs are for positive and negative ip, respec- 
tively, and q = V3/2. Reinserting this result into the free 
energy yields f a (ip). 

The equilibrium densities of the two phases as a func- 
tion of e can then be found by the common tangent con- 
struction, which is equivalent to the requirement that 
the chemical potential /i and the grand potential density 
u> = f — flip must be equal in both phases, 



dip 






fJ-l = Mcq 



W s = fs(lps) - VspJs =U>t = fl{lpl) - Mlpl 



(12) 



(13) 



The solution of these equations yields the equilibrium 
densities as a function of temperature, ?/>f q (e) and ^f q (e). 
The phase diagram of the PFC model exhibits a critical 
point. The parameter e plays the role of an undercool- 
ing, that is, higher e correspond to lower temperatures. 
Furthermore, the phase diagram is symmetric in ip and 
hence exhibits two coexistence zones. We choose for all 
our simulations negative values of ip; for this branch of 
solutions the solid has a higher density than the liquid. 



Finally, for values of ip close to zero, an additional striped 
(nematic) phase can exist, which is not of importance for 
the present work. 

The one-mode ansatz gives a good approximation for 
the phase diagram as long as e remains small. However, 
it turns out that this approximation is not sufficient for 
our purpose since we want to determine excess free ener- 
gies due to surfaces, which requires an excellent precision 
of the bulk values. Therefore, we obtained the function 
f s (ip) from the numerical solution of the free energy min- 
imization for a periodic hexagonal pattern, and used this 
function to perform the common tangent construction, 
which leads to very precise values of [i eq , i/>f q , and tp° q . 

The solid-liquid interfaces have been studied in detail 
in the PFC model and in a Ginzburg-Landau model 33 . 
The main result that is important for the present work is 
that for e small enough, the interfaces are smooth. That 
is, the amplitude of the density waves varies from the 
solid to the liquid over a distance 6 s i that is much larger 
than the spacing between density peaks in the solid. This 
makes it possible to use a multi-scale expansion and to 
obtain a good approximation for the surface tension and 
the order-parameter profile. However, as for the bulk 
densities, this approximation is not precise enough for 
the purpose of the present work. Therefore, the surface 
tension is extracted from the numerical calculations as 
detailed below. The interface thickness 5 s i is obtained 
from a fit of the density profile tp x (y) (the density aver- 
aged over the x direction, which is parallel to the inter- 
face) with a hyperbolic tangent, 



ipx(y) 



^ eq +vr , t 



tr 



tanh 



y 



(14) 



as shown in Fig. [TJ For e = 0.1 (which is used in all 
the simulations in this work), a value of 6 s i « 12.5 is 
obtained. 



Numerical methods 



The standard equation of motion of the PFC model 

.19.20 



d t *P 



6$ 



= (l-e)V^V + 2V> + V b V + V> , (15) 

which reflects the fact that the density field is a lo- 
cally conserved quantity. This equation can be efficiently 
solved by using a semi-implicit pseudospectral formula- 
tion, as detailed in Appendix A. 

However, for the purpose of finding the equilibrium 
states, this is not an efficient method. The reason is that 
the solid and the liquid have different densities, which 
have to be adjusted to their equilibrium values in the 
course of the simulation. Since Eq. (j 1 5[) implies that mass 
is transported by diffusion only, the equilibration time 
scales as the square of the system size. Instead, a more 



rapid numerical scheme can be used, in which tp is treated 
as a locally non-conserved order parameter, while global 
mass conservation is ensured by a Lagrange multiplier. 
The advantage of this "nonlocal" method is that the mass 
can be transported faster since it can be taken at some 
space point and placed at another, as favored by the free 
energy. 

The equation of motion for the nonlocal dynamics is 
derived in Appendix A and can be written as 



d t i> 



8T_ 

Jhp 



= [(e - 1) - 2V 2 - V 4 ] V - ip 3 



fx. 



(16) 



where the Lagrange multiplier li is obtained as 



1 



IX 



L>tX>i, 



[(l-e)^(x)+^ 3 (x)]dx, (17) 



where L x and L y are the side lengths of the rectangular 
simulation box. The Lagrange multiplier is the thermo- 
dynamic chemical potential of the system. In the scheme 
outlined above, the total mass of the system is conserved, 
and the chemical potential evolves with pseudo-time until 
it reaches its stationary equilibrium value. 

Finally, li can also be fixed, and the constraint of global 
mass conservation released. This corresponds to a situa- 
tion described by the grand canonical ensemble, and li is 
the externally imposed chemical potential. The equilib- 
rium state can be reached even faster in this way since 
each point of the system can directly exchange mass with 
the "mass reservoir" . The equation of motion is identi- 
cal to Eq. (|16[) . except that now li is an external pa- 
rameter and independent of time. This method is much 
faster than the others and will be used for almost all of 
the simulations presented below. However, it should be 
emphasized that it is not suitable to simulate isolated 
solid-liquid interfaces or two-phase states within the co- 
existence region, since for such states the density ip is not 
a unique function of li. Such states have therefore to be 
calculated with fixed total mass. 

Here we neglect the effect of thermal fluctuations that 
is traditionally incorporated in the PFC model through 
the addition of a Langevin noise term in the evolution 
equation for the density field, with the amplitude of the 
noise determined by a standard fluctuation-dissipation 
relatio n 19 i 20 ' 38 . This choice is motivated by the fact 
that we focus primarily on computing quantitatively the 
excess interfacial free energies of dry and wet equilib- 
rium grain-boundary states that correspond to stable or 
metastable free-energy minima. This requires an accu- 
rate computation of the free energy of the system that 
is readily obtained from a static crystal density field us- 
ing the "bare" free-energy functional defined by Eq. |(7J), 
but that is considerably more difficult to obtain when 
noise is present. In the latter case, the additional en- 
tropy generated by the fluctuations of the crystal den- 
sity field needs to be computed explicitly to obtain a 
"renormalized" free-energy functional, which is needed 



to compute in a thermodynamically self-consistent way 
the disjoining potential. While such a computation is in 
principle possible (although it would require long simu- 
lations for statistical averaging) it appears unnecessary 
for the computation of static equilibrium properties since 
the bare free-energy functional is derived from a mean- 
field classical density-functional theory framework that 
already contains the effect of microscopic fluctuations on 
the atomic scale. From this fundamental viewpoint, noise 
in the PFC model can only be meaningfully defined in 
the framework of a long-wavelength hydrodynamic the- 
ory where it only acts on length scales larger than the 
correlation length. PFC simulations with noise in this 
hydrodynamic limit and without noise should give essen- 
tially identical results as far as static equilibrium prop- 
erties are concerned. One possible exception is the case 
where different grain-boundary states (corresponding to 
the isolated and paired liquid pool structures already 
mentioned in Sec. lIC[) are separated by small free-energy 
barriers. While such barriers are present in the bare free- 
energy landscape studied here, they could potentially be 
reduced or eliminated in the renormalized landscape due 
to frequent thermally activated transitions between these 
two states. 

The boundary conditions have to be treated with some 
care. The solid phase in the PFC model has a periodic 
structure and can support strain through a variation of 
the wavelength with respect to the equilibrium value. 
However, this variation alters the free-energy density of 
the solid phase. In order to recover the correct equilib- 
rium values in the thermodynamic limit of large system 
size, it is important to ensure that the solid far from the 
grain boundaries is free from strain. Since we use peri- 
odic boundary conditions in both x and y directions, the 
size of the simulation box has to be carefully adjusted 
to contain exactly an integer number of unstrained unit 
cells; this is detailed in Appendix B. 

The initial conditions used to simulate grain bound- 
aries are two solid slabs which are rotated by an an- 
gle = 9/2 in opposite directions. The solid is cre- 
ated using the density field in the one-mode approxima- 
tion ip s (x, y) as given in Eq. (|10p. The solids are initially 
separated by macroscopically large liquid films, where 
ip = ipi- Note that due to the periodic boundary con- 
ditions and the symmetries, there are always two equiv- 
alent grain boundaries in the system. To obtain "dry" 
grain boundaries, ip (or /i in the case of grand canoni- 
cal simulations) is chosen to be within the solid phase. 
Then, in the beginning of the simulations, the liquid 
rapidly solidifies and the grain boundary builds up. Be- 
fore extracting the grain-boundary properties, the system 
is evolved for a much longer time. The approach to equi- 
librium can be monitored by determining the maximum 
difference between the local chemical potential given by 
Eq. ([8]) and the thermodynamic chemical potential (the 
Lagrange multiplier in Eq. (|17[) for conserved total mass, 
or the externally imposed value for grand-canonical sim- 
ulations). 



For = 0, a single crystal is obtained after the liquid 
has disappeared. Due to the symmetry of the hexago- 
nal structure, this happens also when = 30°, but the 
two configurations differ. In the former case, the close- 
packed rows of density peaks are aligned with the x axis 
and hence parallel to the initial liquid layer, whereas in 
the second case, they are aligned with the y direction 
and hence perpendicular to the liquid layer. Therefore, 
configurations with close to or 30° correspond to 
symmetric tilt grain boundaries of inclination <p = 0° 
and <p — 30°, respectively. Furthermore, the misorienta- 
tion is given by 9 = 20 for <p = 0°, but by 9 = 60° - 20 
for <p> — 30° . We recover of course the well-known fact^ 
that there are two equivalent descriptions for each grain 
boundary. In the following, we will investigate the whole 
range of angles < < 30°, which includes low-angle 
grain boundaries of both inclinations. 



III. DETERMINATION OF THE 
GRAIN-BOUNDARY PROPERTIES 

A. General framework 

Experiments and MD simulations are mostly carried 
out at constant temperature, pressure, and total number 
of atoms. Therefore, the appropriate thermodynamic po- 
tential is the Gibbs free energy. In the PFC model, the 
starting point is a Helmholtz free-energy functional. Sim- 
ulations carried out at fixed total mass correspond hence 
to constant temperature (here, e), volume, and particle 
number, and lead to a minimization of the functional 
T . In contrast, if the constraint on the total mass is 
relaxed and the chemical potential is fixed, we have con- 
stant temperature, chemical potential and volume, and 
the relevant thermodynamic potential which is minimized 
by the dynamics is the grand potential, 



Vl = T - n I ip. 
'v 



(18) 



Like the Gibbs free energy, it depends on two intensive 
variables (temperature and chemical potential). We will 
formulate all the subsequent discussion in terms of the 
grand potential, and briefly discuss below how our meth- 
ods and results can be translated to the (N, p, T) ensem- 
ble and the Gibbs free energy. 

The grand potential depends on the intensive variables 
T (here, e) and fi. We will assume in the following de- 
velopments that T is kept constant and that only \x is 
varied. The motivations for this choice will be discussed 
below. Since we have chosen the side of the PFC phase 
diagram where the solid has a higher density than the 
liquid, increasing the chemical potential with respect to 
the coexistence value favors the solid phase. Therefore, 
increasing the chemical potential is analogous to decreas- 
ing the temperature. 

The grain boundary is described as a thin film of liq- 
uid sandwiched between two solids, and the total grand 
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potential of this two-phase system is written as 

Q(/x) = L x [(L y - w)u a (ji) + wwiifi) + 2 7si + V(w)] , 

(19) 
where L x is the length of the grain boundary contained 
in the box (the equivalent of the total surface of grain 
boundary in three dimensions), L y is the system size 
in the direction normal to the grain boundary, w is the 
thickness of the liquid film, and lu s (/i) and LOi(p) are the 
grand potential densities of the bulk solid and liquid, re- 
spectively. Equation (fT9]) is the direct analog in the grand 
canonical ensemble of Eq. ([I]) in the Gibbs ensemble. 

As already described in Sec. I, the last two terms in 
the brackets on the right-hand side describe the excess 
grand potential that is due to the presence of surfaces: 
7si is the surface free energy of an isolated solid-liquid 
interface, and V(w) is the disjoining potential, which de- 
scribes the fact that two solid-liquid interfaces start to 
interact when the distance between them becomes com- 
parable to the range of the interatomic potentials. Since 
V(w) describes the interaction between interfaces, it has 
to tend to zero for well-separated interfaces, V(w) — > 
when w — ► oo. For the form of the disjoining poten- 
tial that has been assumed in the sharp-interface picture 
in Eq. @, a distinction can be made between "attrac- 
tive" grain boundaries for which 7° b — 27 s i < (one 
grain boundary is more favorable than two solid-liquid 
interfaces), and "repulsive" or "wet" grain boundaries 
for which the opposite is true. 

This terminology can be further motivated by defining 
the disjoining pressure II, frequently used in the physics 
of wetting and thin liquid films 18 , 



n = — 



i on 

L x dw 



u>i 



V'(w). 



(20) 



The disjoining pressure has two contributions. The first 
is of thermodynamic origin and changes sign at the melt- 
ing point. Indeed, the grand potential density can be 
expanded in \i around the melting point, which yields 



UJ s -UJlK, - (V>s q - ^r) (M ^ Meq), 



(21) 



where we have used the identity du/dfi = — ip and the 
fact that uj s = oji at coexistence. The second contribution 
in the disjoining pressure arises from the interaction of 
the interfaces. For the simple exponential form of the 
disjoining potential given in Eq. ^, its sign depends 
only on the quantity A7 = 7° b — 2j s \. 

These considerations yield an alternative and quite in- 
tuitive picture of the phenomena already discussed in 
Sec. I. When both contributions of the disjoining pres- 
sure are negative (p, > /i cq , attractive interfaces) the 
film thickness vanishes (w = 0). When both are pos- 
itive (/i < /i cq , repulsive interfaces), the film thickness 
becomes infinite. The more interesting scenarios arise 
when the two contributions arc of opposite signs: for 
attractive interfaces, metastable solids separated by a 
thin liquid film can exist for fx* < \x < /i cq . For repul- 
sive interfaces, finite liquid films exist for /if, > \x > /i eq 



since the repulsion between interfaces competes with the 
thermodynamic force "pushing" the two solids together. 
Here, \x* and Hb are the equivalents of the "breaking" and 
"bridging" temperatures T* and T& defined in Sec. I. 

We would like to point out that the notations used in 
Eq. [2] can easily lead to confusion because of the use of 
the "grain-boundary energy" 7^ in the expression for 
A7. Indeed, the grain-boundary energy of any grain 
boundary, be it "dry" or "wet" , is defined as the total 
excess grand potential per unit length of grain bound- 
ary with respect to a single-phase solid. Therefore, the 
grain-boundary energy is 



Isbip) 



Cl((i) - L x L y uj s (ii) 

Lj; 

(Ul - LO s )w cq (fl) 

+ 2 lsl + V{w cq {fi)), 



(22) 

where the equilibrium film thickness for given chemical 
potential, w eq (ii), is obtained from the condition that 
w eq minimizes the grand potential (which corresponds to 
a vanishing disjoining pressure), 



V'(w cq (n)) = w„(n) - u>i{n). 



(23) 



It can be easily seen that j g b = 7° b only when w = 0. 



It should be emphasized that Eq. (|2"2"]) is completely gen- 
eral, and not limited to the special case of an exponential 
disjoining potential. This relation, which shows that the 
grain-boundary energy and the disjoining potential are 
not independent, can actually be exploited to determine 
the disjoining potential, as will be detailed below. 



B. Liquid film thickness 

To proceed, we need a way to extract the liquid film 
thickness from our simulation data. When the two solid- 
liquid interfaces are well separated, it is easy to define a 
film thickness by the distance between the midpoints of 
the diffuse interfaces. However, this definition becomes 
obsolete when the diffuse interfaces overlap. Another def- 
inition is needed; we choose here to use a Gibbs construc- 
tion. 

When the liquid film is macroscopically large (that is, 
the separation between the two solid-liquid interfaces is 
much larger than the intrinsic interface width), the in- 
terfaces do not interact (the disjoining pressure vanishes) 
and we are in the case of two-phase coexistence, which 
implies that /1 = /i cq (e). The volume fractions of liquid 
and solid are related to the total mass of the system by 
the lever rule. For a one-dimensional system of length 
L y and a film of thickness w, we have 

$L v =<ij) l w + >ip s {L v -w). (24) 

with tyi = $1^ and 'tps = V^ q . 

This is no longer valid when the interfaces interact: 
the disjoining pressure modifies the equilibrium chemi- 
cal potential. However, volume fractions can still be de- 
fined starting from the consideration that the solid is a 



bulk phase which occupies a macroscopic volume. Con- 
sequently, the relation between its density and chemical 
potential is the same as that for a homogeneous solid. 
In contrast, the "liquid" film is microscopic, and hence 
this region does not have the properties of a bulk liquid. 
This is illustrated in Fig. [I] where we show a numerically 
calculated equilibrium state together with a plot of the 
density averaged over the direction parallel to the grain 
boundary. The density exhibits a "dip" and approaches 
the value of the liquid when the film thickness is rela- 
tively large. It exhibits an oscillatory behavior for more 
"dry" grain boundaries, but the average density in the 
grain-boundary region is still different from the one in 
the bulk solid. 

This density change in the grain-boundary region can 
be exploited to define a film thickness. An excess mass 
per unit length of grain boundary can be defined by sub- 
tracting the mass of the homogeneous solid at the same 
chemical potential from the actual mass contained in the 
system, 



^CXC(P) =Ly[lj>- 1p a ((J,)] 



(25) 



Furthermore, it is easy to obtain the density of a bulk 
liquid at the same chemical potential, ipi()J.) from the 
curve of fi{$). Then, the film thickness can be defined 
by the requirement that the density difference of the bulk 
phases times the film thickness is equal to the excess 
mass, 



w [ipi(ji) - ipsifi)] = </w(m)- 



(26) 



Putting these two equations together, we obtain again 
the lever rule, but this time with ip s (fJ>) and tpiif 1 ) m ~ 
stead of the coexistence values. With this definition, the 
film thickness can be extracted with good precision from 
simulations either at fixed total mass (p, is measured in 
the simulation) or at fixed chemical potential (the total 
density ip is measured). 



C. Grain-boundary energy and disjoining potential 

It turns out that the direct numerical determination 
of the grain-boundary energy requires some care. It is 
defined as the excess of grand potential. Contrary to the 
mass excess defined above, the straightforward method of 
subtracting the grand potential of a homogeneous bulk 
solid from the total grand potential of the simulated sys- 
tem leads to large numerical errors. This is most likely 
due to the evaluation of the gradient contributions in the 
free energy. A more precise method is to exploit the de- 
pendence on system size. By dividing Eq. p9[) through 
L x L y and using the definition of the grain boundary en- 
ergy, we obtain that the total grand potential density 
varies with system size at fixed chemical potential as 
u> = uj s + 7 g b/ L y . The grain-boundary energy can there- 
fore be obtained from a plot of lo versus the inverse sys- 
tem size. 
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FIG. 1: (Color online) Density profiles of solid- liquid interface 
and grain boundaries. The complete two-dimensional density 
field ijj(x,y) is shown. The superimposed line gives, for any 
point along the direction normal to the grain boundary, the 
value of the density averaged over the direction parallel to the 
grain boundary (ipx(y) = (1/L X ) J ip(x,y) dx). Top: Solid- 
liquid interface, 6 = 21.8°. Middle: Solid-solid interface close 
to the melting point, tp - -0.1980, 9 - 32. 2°. Bottom: Solid- 
solid interface far from the melting point, tfj = 0.180, 9 = 
32.2°. 
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A second, and slightly simpler way, to obtain the grain 
boundary energy is to perform simulations at a fixed to- 
tal density %j) = i/>o and to use the free energy density, 
which can be directly obtained from the simulations. In- 
deed, since the density in the grain boundary is differ- 
ent from that in the bulk, for a fixed total density and 
length of grain boundary, the bulk density in the solid 
ip s (and therefore also the chemical potential) vary with 
the system size. In the limit L y — > oo, the bulk den- 
sity ip s tends to i/jq and the chemical potential tends to 
the value corresponding to a solid at that density. Us- 
ing the result obtained above, lu = u) s + j g b/L y , and the 
definition u> = f — flip we obtain 



/ = fs(ips) - n{i>a -1P0) + 



Tgfc 
L' 



(27) 



Expanding f s in ijj around ipo, using that df/dip = fj,, 
and inserting the result in the above equation, all the 
terms involving jjl cancel out, and finally we obtain 



f = f.($o) 



Tgb 

Ly 



(28) 



Therefore, we determine the excess grand potential by 
performing simulations at fixed ipo and calculating the 
free energy directly from the free energy functional. Note 
that 7 g b depends on /i and hence also varies with system 
size. However, this gives rise to terms in / that are of 
order 1/Ly and should therefore be small. A plot of the 
total free-energy density versus 1/L y) as shown in Fig. [21 
is indeed well fitted by a straight line. Therefore, we 
can extract the slope and intercept, which correspond to 
7 g b and to the free-energy density in the thermodynamic 
limit, respectively. A numerical error can also be esti- 
mated from the fit if more than two different lengths are 
simulated. The same procedure is also used to determine 
the solid-liquid surface tensions. For this, it is sufficient 
to choose an average density which leads to macroscopi- 
cally large liquid films. 

The disjoining potential can then be obtained in two 
ways. We remark that uii — u> s is a function of \i only, 
and thus V'(fi) is a known function of /i which depends 
only on bulk thermodynamics. The extraction of the liq- 
uid layer thickness from the simulations yields w cq (fi). 
The two can be combined to yield V'(w) which can then 
be integrated to V(w). Alternatively, once the grain- 
boundary energy is calculated, Eq. (|2"2")1 can be used to 
obtain V(/j,), which can again be combined with w(fi) 
to yield V(w). While this second approach avoids a nu- 
merical integration, it is also more costly since for each 
value of the grain-boundary energy several simulations 
with different system sizes have to be performed. 



D. Thermodynamic consistency 

It is useful to comment here on two important points 
with regards to the thermodynamics of interfaces and 
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FIG. 2: Symbols: Free energy densities / of systems with 
the same misorientation and average densities %p but different 
lengths L y perpendicular to the grain boundary, plotted ver- 
sus 1/Ly. Line: Linear fit to the data. The slope gives twice 
the grain-boundary energy, where the factor of 2 is due to 
the fact that in periodic systems there are always two grain 
boundaries. In this example, e = 0.1, ?/> = —0.1, and the mis- 
orientation is 8 — 6°. 



grain boundaries. The first is that with our definition of 
the film thickness, the disjoining potential and the grain- 
boundary energy are entirely thermodynamically consis- 
tent. To show this, let us first remark that Eq. (|2"2")l for 
the grain-boundary energy formally depends on two vari- 
ables, n and w. Taking the differential of this equation, 
we find 



dj gb = [wi—U 
= -Udw 



+ V'(w)]dw + w 



dun 



Qui 
<9/i 



djj, 

(29) 



where we have used the definitions of LT and w and the 
fact that duo/d^i = —ip to obtain the second equality. 
It is clear that the film thickness w plays, for the ex- 
cess free energy, the same role as the volume in bulk 
thermodynamics. Furthermore, since at equilibrium the 
disjoining pressure vanishes, II = 0, the variation of the 
grain-boundary energy is consistent with the fundamen- 
tal definition of interfacial excess quantities**. Indeed, 
since j g b is an excess of grand potential, we can write 



7g6 = /c: 



/#e 



(30) 



where / e xc is the excess free energy; differentiation with 
respect to \x yields 



dn 



-1pe 



(31) 



As a corollary, once a value of 7 g b is known for a sin- 
gle value of fi, the curve 7 g &(^) can be obtained by inte- 
grating the function — ip oxc {n) extracted from the simula- 
tions. Furthermore, formally the grain-boundary energy 
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can also be obtained by keeping /x fixed and integrating 
the "mechanical work" — Hdw over w, noticing that the 
disjoining pressure is non-zero if w 7^ ui oq . This proce- 
dure, however, cannot be carried out in practice since the 
configurations with II ^ are not equilibrium states and 
hence cannot be obtained in simulations. 

The second remark concerns the generalization of our 
definition and procedures to other variables and ensem- 
bles. The various relationships between w, V(w), and 
7 g b obtained above all make use of the fact that the film 
thickness has been defined by a Gibbs construction using 
the interface excess of the density, which is the extensive 
quantity conjugate to the externally controlled intensive 
variable /i. Equivalent constructions can of course be 
performed with other pairs of variables. For instance, in 
their lattice-gas study, Kikuchi and Cahn kept the chem- 
ical potential constant and varied the temperature, while 
they defined the thickness of the liquid layer by the ex- 
cess of entropyii. Similarly, in the (N,p,T) ensemble, a 
film thickness can be defined via the excess entropy for 
varying temperature, or via the excess volume for vary- 
ing pressure. Since, in this ensemble, the volume is no 
longer constant, instead of volume densities as above, 
quantities normalized by the particle number have to be 
used. Nevertheless, following the ideas in Ref.— to treat 
this change in normalization, all the relations given above 
can be translated without difficulties. 

A more complex situation arises if both the chemical 
potential and the temperature (here, e) are allowed to 
vary. For clarity of exposition, we use in the remainder of 
this subsection the temperature T instead of the dimen- 
sionless quantity e. The definition of the grain-boundary 
energy and its variation become 



7gb 
dlgb 



T.Sr 



H4>e 



-Udw 



^exc^M, 



(32) 
(33) 



respectively, where e CX c and s oxc are the interfacial ex- 
cesses of the internal energy and the entropy, respectively. 
Since 7 g f, as well as all the other excess quantities are de- 
fined as excesses with respect to a bulk thermodynamic 
potential that depends on /1 and T, they are all state 
functions, that is, unique functions of the two intensive 
variables fi and T. The same is true of the film thickness 
w, which is defined through an interfacial excess quantity. 
In contrast, the disjoining potential is not a state func- 
tion. This is easy to see when considering the equilibrium 
condition for the film thickness, Eq. (|23|) : its right-hand 
side, u s — u>i, now depends on the two variables \x and 
T, which implies that V'(w) has the same dependency. 
If this is to be integrated to a function of a single vari- 
able w, a direction in the space spanned by fi and T 
has to be specified. Another way to state the same fact 
is to remark that in Eq. (|22p , j g b depends on the two 
independent variables /i and T, whereas the "reference 
value" 2^/ s i depends only on one independent variable, 
since 7^ is only defined on the coexistence line in the 
phase diagram, /x cq (T). For a given point away from this 
line, where 7 g b is still defined, V can be defined only 



if a reference point on the coexistence line is specified. 
This amounts to specifying the path in the state space 
that is to be followed. In the developments above, we 
have supposed a particularly simple path, namely, a con- 
stant value for one of the variables. It would be possible 
to extract from our PFC model disjoining potentials at 
constant pressure or at constant density: for both cases, 
the bulk equation of state for the solid (which can be 
obtained from / s (/i,T)) fixes a relation between /i and 
T, and V'(w) can be integrated along this path. Note, 
however, that this procedure requires the calculation of 
both the excess mass and the excess entropy. 

In summary, the definition of the disjoining potential 
is only meaningful if the corresponding path in thermo- 
dynamic state space is specified, and the knowledge of 
a single disjoining potential yields only a partial knowl- 
edge about the premelting transition. The more general 
quantity is the grain-boundary energy, which is the ther- 
modynamic potential for the interfacial excess quantities. 
If its dependence with respect to the two intensive vari- 
ables is known, all the possible disjoining potentials can 
be easily extracted using Eq. 



E. Choice of simulation parameters 



Having the discussion of Sec. IIIIDI in mind, we need 
to choose a particular path in the state space to inves- 
tigate the disjoining potential. It would be possible to 
approach the melting transition from the solid side for 
a fixed chemical potential by decreasing e. However, ex- 
tracting the excess entropy is far more delicate than ex- 
tracting the excess mass. Therefore, in the following we 
prefer to keep e fixed to 0.1 (a value that has been ob- 
tained for the equilibrium solid-liquid interfaces in pure 
iron with body-centered-cubic crystal orderin g 33 ! 34 ) and 
to explore the melting transition by varying the chemical 
potential /1. 

For the subsequent presentation of the results, we will 
use the following rescaled variables: 



A = 



1>-^ 



.cq 



V's 



vr 



(34) 



This corresponds to a supersaturation. Furthermore, we 
define a scaled chemical potential by 



Mcq - A* 



dip s 



w 



cq 



^D 



(35) 



where the sign is chosen to stress the analogy between u 
and a temperature^.: for u < (u > 0), the solid (liquid) 
is the favored phase. For \x close to the coexistence value, 
the numerator can be expanded in rp, which yields u ~ 
—A. We list in Table |T] the values of all the quantities 
needed for this scaling. Furthermore, we will often rescale 
the grain-boundary energy by 2^, and lengths by a, the 
lattice spacing of the hexagonal crystal. The values of 
these quantities are also given in Table HI 
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TABLE I: Numerical values of various quantities needed to 
scale the density, chemical potential, grain-boundary energy, 
and lengths, for e = 0.1. 

Quantity Symbol Value 



Solid density at coexistence ?/>s q 

Liquid density at coexistence i/>, eq 

Chemical potential at coexistence /i eq 
Slope of the curve fi versus ip s dfi/d'ips 
Solid- liquid surface tension (x2) 27 s i 
Lattice constant a 



4- 



0.19696406 
-0.2068060 
0.19497015 
0.731218 
0.00192 
7.2552 



IV. RESULTS 
A. Structure of the grain boundaries 

Our simulations reveal that there is a strong differ- 
ence in behavior between high-angle and low-angle grain 
boundaries. In order to illustrate first a few important 
features, we show in Figs. [3] and H] snapshot pictures of 
a high-angle and a low-angle grain boundary of inclina- 
tion <j> — 0°, for different values of u. Furthermore, we 
plot in Fig. [5] the curves of film thickness w versus scaled 
chemical potential u corresponding to the same two grain 
boundaries. We have checked that canonical and grand- 
canonical simulations (fixed total mass and fixed chemi- 
cal potential, respectively) yield identical results for the 
film thickness and the grain-boundary structure. 

For both high-angle and low-angle grain boundaries, 
the film thickness becomes negative far below the melt- 
ing point. Indeed, formally, since the film thickness is 
defined via an excess mass, it does not need to remain 
positive. A negative film thickness corresponds to an ac- 
cumulation of mass in the grain boundary instead of the 
depletion observed in Fig. [TJ When u is increased, the 
film thickness becomes positive, but remains small until 
the vicinity of the melting point is reached. For the high- 
angle grain boundary, the film thickness then increases 
rapidly and diverges as the melting point is approached 
from below; this is the behavior expected for a repulsive 
grain boundary. In the snapshot pictures, it can be seen 
that the liquid film is rather homogeneous, that is, it has 
approximately the same width at every point. 

The low-angle grain boundary depicted in Fig. [3] con- 
sists of individual dislocations separated by distances 
that are larger than a few lattice spacings. Here, the 
"liquid" first appears in the form of "pools" around the 
dislocations, and there is no homogeneous film of liquid. 
Furthermore, as the melting point is approached, a struc- 
tural transition occurs: the dislocations form pairs; that 
is, two dislocations join and are surrounded by a common 
liquid pool. This transition is accompanied by a jump in 
the film thickness w. Furthermore, this structure can be 
"overheated"; that is, such states exist even for u > 0, 
which indicates an attractive grain boundary. The pools 
grow in size, thus reducing the strength of the "bridges" 
of solid. At a critical overheating, the bridges break, and 
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FIG. 3: (Color online) Snapshots of a high-angle grain bound- 
ary with 6 — 32.2° for different values of u, which increase 
from bottom right to top left (see text and Fig.[5]for details). 
The "liquid" forms a rather homogeneous film. Only part of 
the simulation box is shown. 



the whole system becomes liquid. 

Most of the features described above - dependence of 
the film thickness on u, transition from attractive low- 
angle to repulsive high-angle grain boundaries, existence 
of overheated states - are also present for the symmetric 
tilt grain boundaries of inclination <j) = 30°. However, for 
this inclination there is no transition from single disloca- 
tions to dislocation pairs. This transition was also not 
observed in the three-dimensional PFC study with bcc 
symmetry in Ref. 26 . It can hence be concluded that its 
occurrence depends on the detailed microscopic structure 
of the grain boundary. 

Let us now give a more detailed description of the tran- 
sition between the high-angle and the low-angle regimes. 
In Fig. [S] we show the curves of w versus u, for various 
misorientations, in the vicinity of the melting point, for 
the two inclinations = 0° and <f> = 30°, respectively. 
We recall (see the discussion in Sec. Ill B| ) that due to the 
hexagonal symmetry the two curves shown for = 0°, 
9 = 32.2° and for <\> = 30°, $ = 27.8° actually describe the 
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FIG. 4: (Color online) Snapshots of a low-angle grain bound- 
ary with 6 = 9.4° for different values of u, which increases 
from bottom right to top left (see text and Fig.[5]for details). 
The grain boundary consists of individual dislocations and 
undergoes a structural transition. Only part of the simula- 
tion box is shown. 
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FIG. 5: Ratio of film thickness w to lattice spacing a as a 
function of scaled chemical potential u for two different grain 
boundaries. The inset shows a blowup of the vicinity of the 
melting point. The symbols mark the states that are depicted 
in the snapshot pictures in Figs. [3] and [4] 



FIG. 6: (Color online) Film thickness raasa function of u for 
various misorientations, for symmetric tilt grain boundaries 
of inclination 4> — 0° (top) and 4> — 30° (bottom) close to 
the melting point. All angles are given in degrees, and a 
vertical line has been drawn at the melting point u = 0. Inset: 
Film thickness of the three largest angles versus — ln( — u): the 
divergence of the film thickness is logarithmic. 



same grain boundary. All curves have been calculated by 
simulations at fixed chemical potential. The final state 
of a given run was used as initial condition for the next 
one at slightly different chemical potential. 

The insets of Fig. [6] show the film thickness for the 
three largest misorientations for both inclinations, all cor- 
responding to repulsive interfaces, versus — ln(— u). For 
large film thickness, the curves become linear, which is 
the dependence that is expected for an exponential dis- 
joining potential from Eqs. @ and ([3]). According to 
Eq. (J3J), the slope of this linear part is the decay length 
8. We find a value of 8 « 5.8, which is approximately half 
of the thickness of the solid- liquid interfaces 8 s i ss 12.5, 
and comparable to the wavelength of the dominant den- 
sity waves of the hexagonal structure (which is equal to 
2ir in our scaling). 
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It can be seen that the transition between repulsive 
and attractive behaviors occurs at an angle of 9 C rs 14° 
for both inclinations. This transition is smooth in the 
sense that the critical value u* where the solid bridges 
break decreases with increasing misorientation and seems 
to tend to zero at the transition angle without exhibiting 
a jump. Furthermore, the thickness of the liquid layer 
at the melting point, w m — u>(0), increases with misori- 
entation and seems to diverge continuously when 9 C is 
approached from below. The precise nature of this diver- 
gence remains undetermined. Its detailed study would 
require simulations in a narrow range of misoricntations 
close to the critical angle, which is quite cumbersome be- 
cause of the geometrical constraints that arise from the 
finite size of the simulation box. 

It is important to stress that this transition does not 
coincide with a structural transition of the grain bound- 
ary. The curves of w versus u for 9 = 32.2° and 9 = 17.9° 
in Fig. [6] are very similar; however, the structure of these 
grain boundaries is quite different. In all the snapshot 
pictures in Fig. [3j the grain boundary is a plane of mir- 
ror symmetry for the density field. This is not the case 
for 9 = 17.9°: far from the melting point, this grain 
boundary consists of individual dislocations such as the 
low-angle grain boundary shown in Fig. 0J The tran- 
sition from single dislocations to dislocation pairs also 
occurs, but far from the melting point, around u = — 1.4. 
When the melting point is approached, a continuous tran- 
sition from a state similar to the uppermost left picture 
in Fig.|3]to one that looks like the uppermost left picture 
in Fig. [3] occurs: the liquid pools around the dislocation 
pairs increase in size and finally merge to give rise to a 
fairly homogeneous film. The "liquid pools" separated 
by "solid bridges" are therefore present in the vicinity of 
the melting point both for repulsive and attractive grain 
boundaries. 

In Fig. [SI jumps in the film thickness can be seen in 
the curves for 9 = 8.6°, 9.4°, and 9 = 10.3°; they cor- 
respond to the occurrence of the transition from single 
dislocations to dislocation pairs. The value of u at which 
this transition occurs increases with decreasing misori- 
entation. As mentioned before, for 9 — 17.9°, it occurs 
far below the melting point; for 9 = 8.6°, it occurs only 
above the melting point. For an even lower misorienta- 
tion, 9 — 6.0°, it does not occur at all: the liquid pools 
around the single dislocations increase in size until the 
solid "bridges" between them break. 

This structural transition exhibits a hysteresis. In 
Fig. [7] we show two curves of w versus u for 9 = 9.6° 
that are computed in different ways. In the first, we start 
from a single-dislocation state at low values of u and then 
perform successive simulations with increasing u, taking 
the final state of the previous simulation as initial con- 
dition. In the second, we start from a dislocation-pair 
state and successively decrease the value of u. It can be 
seen that there exists a range of u in which both single- 
dislocation and dislocation-pair states are stable. This 
indicates that there exist at least two distinct branches 
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FIG. 7: Various grain-boundary states for inclination c/> = 0° 
and misorientation 6 = 9.4°. The curve w versus u ex- 
hibits a hysteresis: filled symbols are calculated with in- 
creasing u (same curve as shown in Fig. |6}, and open sym- 
bols with decreasing u. The snapshot pictures show differ- 
ent grain-boundary states, all at U — —0.072. Top: single- 
dislocation state, middle: dislocation pair calculated with a 
single-dislocation state as initial condition, bottom: disloca- 
tion pair calculated with the initial conditions described in 
SecHTBl 



of grain-boundary solutions. In addition, there exist dif- 
ferent configurations for the dislocation-pair state, as is 
shown in the snapshot pictures in Fig. if the simu- 
lation is started from the initial conditions described in 
Sec. lIIEfl the dislocation pair state exhibits a mirror sym- 
metry with respect to the plane of the grain boundary, 
as also seen in the snapshots in Fig. 2J The dislocation- 
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pair states obtained starting from a single-dislocation 
state do not exhibit this symmetry. However, the size 
and shape of the "liquid pool" are quite similar, and the 
film thickness extracted from the two different configu- 
rations is identical up to the numerical precision. We 
have also found that in the case of the low-angle grain 
boundary that does not exhibit the transition (0 = 6.0°), 
dislocation-pair states can be obtained starting from the 
initial condition in Sec. Ill Bt they form a second branch 
of solutions for this misorientation. 

The hysteretic nature of this transition and the depen- 
dence of the final state on the initial conditions are clear 
indications that the free-energy functional of the PFC 
model has several distinct minima which correspond to 
grain-boundary states of different grain-boundary ener- 
gies. For low values of u, the single-dislocation states 
have a lower energy and the dislocation-pair states cor- 
respond to a metastable minimum, whereas the inverse 
is true for high values of u. Since the extraction of 
the grain-boundary energy is delicate, we have not pin- 
pointed the exact value of u where the two states have 
equal energies. But from the general phenomenology of 
hysteretic transitions, it can be expected to lie approxi- 
mately in the middle of the bistable range. Furthermore, 
since our simulations do not include thermal fluctuations, 
the termination of the metastable branches corresponds 
to the disappearance of the local metastable minimum. 

The existence of metastable states raises the question 
of whether other grain-boundary configurations, distinct 
from the ones depicted in Figs.[3]to[71 might exist. We in- 
vestigated this question by performing several runs with 
different initial conditions for numerous parameter sets. 
We did not find any new grain-boundary states in the 
vicinity of the melting point. For states far from the melt- 
ing point, we have occasionally observed distinct configu- 
rations that exhibit differences in the local arrangements 
of the "atoms" around the dislocations and different total 
number of "atoms", which is possible since, even at fixed 
total mass, the total number of "atoms" is not fixed in 
the PFC model. No further investigation of these multi- 
ple states was carried out. 

The curves of w versus u for intermediate misorien- 
tations that consist of dislocation pairs exhibit a verti- 
cal slope at the "break point" . It is tempting to believe 
that the solution branch continues beyond that point and 
bends back to reach u = when w — > oo. This would be 
expected if the state of the grain boundary can be prop- 
erly described by the single variable w; such solutions 
have been found in phase-field studies of grain-boundary 
premelting2£i2£ . In our PFC model, these solutions can- 
not be obtained in simulations at constant chemical po- 
tential, since they are unstable. We have tried to ob- 
tain these states by simulations with fixed total mass. In 
this case, mass conservation yields a constraint on the 
film thickness which should stabilize these states. How- 
ever, our attempts were not successful due to the occur- 
rence of a new instability. Since there are always two 
distinct grain boundaries in our system due to the peri- 



odic boundary conditions, a symmetry breaking can oc- 
cur which leads to the formation of a "thick" and a "thin" 
liquid film instead of two liquid films of equal thickness. 
This indeed happens when the film thickness is larger 
than the value corresponding to the "turning point" . A 
simple explanation for this instability will be given be- 
low. In contrast, all the curves w(u) for low-angle grain 
boundaries with inclination <fi — 30° as well as the curve 
for the lowest misorientation for = 0° (9 — 6.0°) do 
not exhibit a turning point, but break off with a finite 
slope. These grain-boundary states all consist of single 
dislocations. 

From these results it can be concluded that the mech- 
anisms that lead to the breaking of the "solid bridges" 
and to the instability of the overheated solution branches 
depend on the detailed structure of the grain bound- 
ary. Qualitatively, the difference in behavior can be un- 
derstood from geometric considerations. As mentioned 
above, a vertical slope at the break point would be ex- 
pected for homogeneous liquid films that can be faith- 
fully described by a single variable, the film thickness w. 
The elongated liquid pools surrounding the dislocation 
pairs are more similar to a homogeneous liquid film than 
the round liquid pools surrounding single dislocations. 
Thus, it is not surprising that the behavior of the former 
is closer to the one of a homogeneous film. 

It is clear from the above results that the description 
of a grain boundary by a single variable (the thickness) 
is very crude. However, we have found no other obvi- 
ous quantity that could play the role of a supplementary 
state variable. Therefore, for all the following develop- 
ments we will restrict our level of description to the sin- 
gle variable w, leaving a more detailed investigation as 
a subject for further study. Also note that, in princi- 
ple, a distinct grain-boundary energy and disjoining po- 
tential are associated with each of the distinct solution 
branches. To simplify the picture, we will ignore this fact 
and display in the following unique curves for the grain- 
boundary energy and the disjoining potential. Since the 
film thickness w exhibits a jump, V(w) has a "step", 
and the grain-boundary energy a discontinuity in slope. 
These features are, however, so small that they can be 
hardly distinguished in the following plots. 



B. Grain-boundary energy 

We calculate the grain-boundary energy as described 
in Sec. IIII CI by performing simulations at fixed density 
■0 for several different system sizes L y and using Eq. [)28p 
to extract j s t,. In order to keep the presentation consis- 
tent, we will nevertheless discuss the results as a func- 
tion of u, which can be easily obtained for given density 
using the curve /x s (V>). All the data shown in this sub- 
section are for grain boundaries of inclination cj> = 0°. 
The grain-boundary energy is plotted versus misorienta- 
tion in Fig.[5]for various values of the chemical potential. 
Two clear tendencies can be seen. First, for any fixed 
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misorientation, 7 g f, increases monotonously when u de- 
creases. Second, for a fixed supersaturation, 7 g (, increases 
monotonously with the misorientation for small angles. 

The latter dependency can be well understood in terms 
of the Read-Shockley law^S, 



7g6 



Go 



4na(l — a) 



(9 [1 - ln(27T0) + ln(ao/r )] , (36) 



where r$ is the core radius of the dislocations, a is the 
lattice constant, a = \/3/2 is the ratio of the distance 
between close-packed planes and the lattice spacing, G 
is the shear modulus, and a is Poisson's ratio. The 
elastic properties of the PFC model can be determined 
analytically in the one-mode approximatio n 19 ! 20 . The 
resulting elastic constants are Cn/3 = C\i = C44 = 



13-0 — yl5e — 36V' 2 ) /75. The bulk modulus can then 

be calculated to be Y = 2C44 and the shear modulus 
G = C44. Furthermore, the three-dimensional Poisson's 
ratio is a = (3Y - 2G)/[2(3Y + 2G)] = 1/4. An esti- 
mation for the core radius, r — aexp(— 0.5) w 4.4 has 
also been giver* 2 ^. Note that we have used here the stan- 
dard (three-dimensional) version of the Read-Shockley 
law and the elastic constants. Tt can be shown^i that 
this is identical to the two-dimensional expressions given 
by Elder and Grants. Furthermore, this formula differs 
from the standard one for cubic materials by the pres- 
ence of the factors a, which are due to the fact that the 
average distance d between dislocations is d ~ aa/ sin 
(instead of d ~ a/ sin 9 for a cubic material). It should 
also be emphasized that this expression is only valid for 
grain boundaries of inclination = 0°; in the general 
case, several additional terms depending on the inclina- 
tion angle have to be included^. 

For each supersaturation, we fixed the shear modulus 
to its density-dependent analytical value, and performed 
a least-squares fit of our data to the Read-Shockley law, 
with tq as the only fit parameter. Since Eq. (|36|) is only 
valid for small misorientations, for the fit only systems 
with 6 < 15° have been included. It can be seen that the 
fit is excellent. In the inset, the core radius r$ obtained 
from the fit is shown as a function of u. It is almost 
constant and close to the theoretically estimated value 
for large values of |u|, and increases when the coexistence 
region is approached (u — ► 0). 

It turns out that the variation of the grain-boundary 
energy with u will be crucial for the further discussions. 
We recall that this variation is directly linked to the liq- 
uid film thickness by Eqs. (J3T|) and (j26|) . The grain- 
boundary energy is shown as a function of u for three se- 
lected misorientations in Fig. [9] The symbols are values 
that have been directly obtained from simulations with 
varying system size. The full lines are obtained by inte- 
grating Eq. (|3ip , where the integration was started from 
the data point at u = —5.674. It is clear that the two pro- 
cedures give fully consistent results. In the inset, the data 
for 7 g f, obtained by integration are shown in the vicinity 
of the melting point. It should be mentioned that direct 
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FIG. 8: (Color online) Ratio of grain-boundary energy to 
twice the solid-liquid free energy for grain boundaries of incli- 
nation cj> — 0° as a function of the misorientation 9, shown for 
different chemical potentials u. The lowest curve, u = —0.006, 
is very close to the melting point, while the upper curve is far 
inside the solid region. The lines are fits to the Read-Shockley 
law, Eq. (|36p , using the values of j g t for 6 < 15° . The shear 
modulus has been fixed to the theoretical value at the cor- 
responding chemical potential, and the only fit parameter is 
the dislocation core radius 7*0 • Inset: The core radius ro ob- 
tained from the fits as a function of chemical potential. For 
comparison, the lattice constant is a « 7.255 and the value 
estimated by Elder and Gran1)2£ is ro « 4.4. 



calculations of the grain-boundary energy in this regime 
are quite difficult, since the grand potential differences 
between solid and liquid (and hence the driving forces) 
are small, so that long equilibration times are needed. 
The different behaviors of repulsive and attractive grain 
boundaries can be clearly seen. For the repulsive grain 
boundary (upper curve), 7 g b tends to 27 s ; from above. 
Since (/; exc ~ ln(— u) close to the melting point, we have 
7 g 6 — 27 s i ~ u ln(— u) + u, as expected from Eq. §5§ of the 
sharp-interface theory. The curve has an infinite slope at 
u = (corresponding to a diverging film thickness), but 
the logarithmic divergence is too slow to be clearly distin- 
guished in the figure. For the attractive grain boundaries 
(the lower two curves), "f g b becomes lower than 27 s ; be- 
fore the melting point is reached. It continues to decrease 
beyond the melting point, until the metastable solution 
branch ends. 

As mentioned above, the dislocation core radius ex- 
tracted from the fits to the Read-Shockley law is almost 
constant over a wide range of u. An interesting corollary 
of this finding is that the variation of the grain-boundary 
energy with chemical potential can be accounted for al- 
most entirely by the change in the elastic constants. To 
illustrate this point, we show in Fig. [H] as dash-dotted 
lines the predictions of the Read-Shockley law with a 
constant value for the core radius ro = 3.75 for the two 
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FIG. 9: (Color online) Symbols: grain-boundary energy de- 
termined from direct simulations, lines: grain-boundary en- 
ergy calculated by thermodynamic integration, dash-dotted 
lines: prediction of the Read-Shockley law with ro = 3.75 
and elastic constants depending on u. See text for details. 



FIG. 10: (Color online) Disjoining potential for inclination 
4> = 0° and various misorientations (inset: detailed view of 
the region around the minimum). The angles are given in 
degrees. 



lowest misorientations. Clearly, the main variation of 7 g (, 
with supersaturation is well reproduced. The ratio of the 
predicted to the numerical value remains close to unity 
up to u ~ —2 and then increases sharply to about 1.4 at 
the melting point. This is natural since the assumption of 
constant core radius breaks down. The highest misorien- 
tation shown in Fig. [5]is too large for the Read-Shockley 
law to be applicable. But from the figure it is clear that 
the variation of 7 g f> with the chemical potential is very 
similar to the one of the low- angle grain boundaries. It is 
therefore reasonable to assume that this variation is also 
mainly controlled by the elastic constants. 



C. Disjoining potential 



As outlined in Sec. IIIIC1 two methods can be used 
to extract the disjoining potential from the simulation 
data: V'(w) can be integrated using the data of w(fi), or 
V(w) can be directly deduced from the grain-boundary 
energy using Eq. (|22p . Both methods yield consistent 
results that are shown in Fig. [TUJ For the high-angle 
grain boundaries, V(w) decreases monotonously; it can 
be actually quite well described by an exponential func- 
tion as in Eq. @. In contrast, for the low- angle grain 
boundaries, the disjoining potential is non-monotonous: 
starting from a positive value at w = 0, it decreases, falls 
below zero and exhibits a minimum for some intermedi- 
ate values of w. It then starts to increase until it reaches 
the point where the curve w(/j,) terminates. 



V. DISCUSSION 

It is instructive to discuss some aspects of our above re- 
sults in more detail and to compare them with the predic- 
tions of the sharp-interface theory. Three questions are 
of particular interest: what is the interpretation of the 
non-monotonous disjoining potentials, what determines 
the critical angle for the attractive-to-repulsive transi- 
tion, and what can be said about the overheated grain 
boundaries and about transitions between different grain- 
boundary states ? 

The shape of the disjoining potential for low-angle 
grain boundaries can obviously not be described by the 
simple exponential form of Eq. ^. This potential cor- 
responds to a short-range repulsion, but a long-range at- 
traction. We have not found a simple analytical formula 
for this potential; however, a few of its properties can be 
readily understood. For instance, Eq. (|20)l tells us that 
V'(w) = implies lu s — lji, which is only the case for 
jj, = fi cq (u = 0). The minimum of the curve V(w) cor- 
responds therefore to the intersection of the curves w(u) 
with the u = axis in Fig. [6] Since Eq. (|22|) yields, for 
uj s = u>i, V(w) = 7 g b — 27 s ;, this implies that the depth of 
the potential well is given by the difference of the grain- 
boundary energy at the melting point and twice the solid- 
liquid free energy. Furthermore, the value V(0) corre- 
sponds to the grain-boundary energy of a completely dry 
grain boundary 7^. The height of the "repulsive part" 
of the disjoining potential is therefore given by the differ- 
ence between this value and the grain-boundary energy 
at the melting point. Any system in which the grain- 
boundary energy increases with decreasing homologous 
temperature will therefore exhibit a repulsive part in the 
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disjoining potential, even if the grain-boundary is attrac- 
tive at the melting point. In addition, from Eq. (|31[) and 
Eq. (|23|) it is easily seen that the variation of the grain- 
boundary energy with u is proportional to the negative 
of the liquid film thickness. Therefore, the disjoining po- 
tential is repulsive below the melting point for any grain 
boundary that exhibits a finite film thickness. 

It is also easy to show that the points where the curve 
w(/j.) exhibits a vertical tangent correspond to an inflec- 
tion point in the potential V(w). For this, it is sufficient 
to take the derivative with respect to w of Eq. (f2"3"| . which 
yields 



<9/i 9/j, J dw 
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(37) 



The sign of the second derivative of the disjoining poten- 
tial is hence determined by the derivative dfi/dw, which 
is zero at the turning point of w(/j>). As a consequence, 
V(w) is concave for large values of w. This yields a simple 
explanation for the instability that leads to a symmetry 
breaking between the two grain boundaries in the sim- 
ulation box, which was described in Sec. IIV Al at fixed 
density, the sum of the two film thicknesses w% and w% is 
approximately fixed by the lever rule. If w± — u>2 = w, 
and w is located in the concave part of the potential, 
the system can lower its total energy by making one film 
wider and the other one thinner. 

Let us now consider the transition from attractive to 
repulsive grain boundaries. As discussed in Sec. I, the 
sharp-interface theory predicts this transition to occur 
when 7 g t, = 2^^. Now consider the different curves of 
7gb versus misorientation shown in Fig. [5] For the two 
lowest values of u, this curve intersects the line corre- 
sponding to twice 7 s z for a misorientation of 9 « 2°, 
much smaller than the transition angle obtained from 
the curves w(u) in Fig. [SJ However, with increasing u, 
this intersection point moves toward larger angles, and 
for the highest value of u investigated, the intersection is 
at about 9°. Furthermore, it was shown above that when 
the disjoining potential exhibits a minimum, its value at 
this minimum is equal to 7 g {, — 2^ s i at the melting point. 
Since the transition to repulsive grain boundaries occurs 
when this minimum disappears, it is to be expected that 
the correct transition angle is obtained when the crite- 
rion 7gf, = 27 s i is used with the grain-boundary energy 
calculated exactly at the melting point. 

In order to obtain more precise information on this 
question, it would be desirable to have accurate values 
for the grain-boundary energy at the melting point as 
a function of misorientation. However, as already men- 
tioned, it is numerically very difficult to obtain values 
for the grain-boundary energy close to the melting point, 
especially for angles close to the repulsive-to-attractive 
transition, since the grand potential differences between 
the different states become extremely small. The best 
way to obtain reliable data close to the transition is to 
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FIG. 11: Ratio of the grain-boundary energy and the melting 
point 7™ and twice the solid-liquid free energy as a function 
of misorientation for inclination cf> — 0° . The line is a fit to 
the Read-Shockley law for the points with 6 < 10° . 



integrate Eq. (|31|) up to the melting point. The result is 
shown in Fig. 1111 we estimate the error bars for these data 
to be on the order of the size of the symbols. For low- 
angle grain boundaries, the dependency of 7 g & on 9 can 
still be described by the Read-Shockley law. However, 
for higher angles, when the grain boundaries consist of 
dislocation pairs, the dependence of 7 g b on 9 is extremely 
weak: for 9 = 10.4°, which is the first point that clearly 
deviates from the Read-Shockley law, j g b/(2-f s i) « 0.99, 
so that the variation of j g b between this misorientation 
and the first repulsive grain boundary at 14.1° is only 
about 1%. Clearly, it is very difficult to describe pre- 
cisely this regime. In addition, it is not clear whether 
it is generic, since the grain boundaries of inclination 
<f> = 30° do not exhibit the structural transition to dislo- 
cation pairs. 

A very interesting point is that the Read-Shockley law 
is still valid for low-angle grain boundaries, even at the 
melting point. This can be used to obtain a reasonable 
estimate for the critical angle as the solution of the equa- 
tion 



Ga 



47ra(l — a) 



9 C [1 - ln(27r0 c ) + ln(aa/r )] = 2 7s/ . (38) 



However, it is crucial to take into account the variation 
of the grain-boundary energy with chemical potential (or 
temperature). Indeed, the values for the grain-boundary 
energy are usually determined in experiments or atom- 
istic simulations for temperatures far below the melting 
point. As pointed out above, if these values are used 
to predict the critical angle, a completely wrong result 
is obtained. The variation of the grain-boundary energy 
with chemical potential (or temperature) arises from two 
distinct effects: the variation of the shear modulus and 
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FIG. 12: Critical value u* that corresponds to the limit of su- 
perheated states associated with the breaking of solid bridges. 
Symbols: simulation results, line: sharp-interface prediction 
according to Eq. (|40)l . 



the premelting around dislocations, which leads to an in- 
crease of the core radius in the Read-Shockley law as 
shown in Fig. [8] If the "low-temperature" values for 
both G and rn are used in Eq. ([55)1 . we find 9 C rs 2° 
(see Fig. [8]), clearly too low. If only the variation of G 
is included (that is, Eq. (f3"B")) is used with the value of 
the shear modulus at the melting point, but with the 
low-temperature value rn = 3.5 for the core radius), the 
prediction becomes 9 C ss 6°. Finally, the improved es- 
timate corresponding to Eq. ([B]) is obtained when the 
values at the melting point of both the shear modulus 
and the core radius are used, which yields the prediction 
9 C rs 10°. The curve of 7^; that is obtained with these 
parameters is shown in Fig. 1111 Clearly, the limited ac- 
curacy of this final estimate is due to the fact that the 
Read-Shockley law applies to grain-boundary states with 
individual dislocations. Thus, it does not take take into 
account the complex structural changes in the bound- 
ary, which will tend to reduce the grain-boundary energy 
further from its value estimated from the Read-Shockley 
law. Therefore, the value of 9 C obtained from this law is 
likely to be a lower bound estimate of the actual value. 

Let us now come to the estimation of the critical value 
u* that limits the range of "overheated" states. The 
transposition of the sharp-interface prediction, Eq. (0|, 
to our variables is 



u>i(u*) - lo s (u*) = 



A7 
'~5~' 



(39) 



Expanding the grand potential around the melting point, 
and using the definition of u, Eq. (|35[) . we find 



2 7 sz 



da 

dips 



($ 



C q 



fs 



A7 



(40) 



The ratio in brackets on the right-hand side can be calcu- 
lated using the values from TableQ]and the values S « 5.8 
and 27s/ « 0.00192 extracted from our simulations (see 
Fig. [S]). From the preceding discussion, it is clear that 
a reasonable estimate can only be obtained if A7 is cal- 
culated with the grain-boundary energy at the melting 
point. In Fig. [TSJ we plot the values of u* obtained from 
our simulations (that is, the values of u where the curves 
of w{u) terminate) versus — A r )/{2 n j B {), using the values 
for 7 g (, in Fig. QTJ together with the theoretical predic- 
tion. It can be seen that this prediction gives reasonable 
values for the grain boundaries close to the transition 
that consist of dislocation pairs, even if the simulation 
data cannot be well described by a straight line. In con- 
trast, the sharp-interface theory strongly overestimates 
the value of u* for low-angle grain boundaries consisting 
of single dislocations. 

The failure of the sharp-interface theory to predict the 
superheated range of grain boundaries is not surprising 
since the liquid phase domains consist of liquid pools in- 
stead of a thin liquid film of constant thickness as as- 
sumed in this theory. Recently, Berry et al££- developed 
a simple theory of grain-boundary wetting tailored to the 
liquid-pool geometry, which assumes that wetting occurs 
when pools coalesce, or, equivalently, when their radius r 
is equal to half of the distance d between dislocations. By 
calculating the shift of the melting point due to the dis- 
location elastic strain energy, they also obtained the scal- 
ing relation u* ~ — (a/cf) 2 , where the dimensionless pro- 
portionality constant is related to the clastic constants. 
This scaling relation, together with the geometrical coa- 
lescence condition r = d/2 — a/ (2 sin 9), yields the pre- 
diction u* ~ — sin 2 9. As shown in Fig. [T3J our results 
for the grain boundaries of the <f> = 30° inclination that 
consist of unpaired liquid pools indeed show that u* is 
reduced by an amount proportional to sin 9, consistent 
with this prediction. One important difference, however, 
is that the 9 — intercept of the curve u* (9) is finite in 
our simulations, consistent with the existence of super- 
heated metastable grain boundaries, while the theory of 
Berry et al~& predicts that liquid pools always coalesce 
below the melting point (u*(0) = 0). We expect u*{9) 
to be generally positive in the limit of vanishing misori- 
cntation since a finite bulk thermodynamic driving force 
favoring the liquid phase is necessary to overcome the 
nucleation barrier imposed by the solid-liquid interfacial 
energy, which remains finite even in the presence of elas- 
tic strain energy around the dislocation cores. 

It is interesting to note that the linear interpolation of 
the data in Fig. [T3] predicts that u* vanishes for sin 2 9 « 
0.06, which corresponds to a misorientation of 14.6°, in 
good agreement with our previous estimate for 9 C . For 
misorientations above this value, no overheated states 
can exist, and there is hence no discontinuous transition 
between "dry" and "wet" grain-boundary states. 

To further test this theoretical picture for the <f) = 30° 
inclination, we have extracted the pool radius from the 
data for the film thickness using the simple geometrical 
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FIG. 13: Critical value u* versus sin 2 9 for low-angle grain 
boundaries of inclination cf> — 30°. Inset: Ratio of the liquid- 
pool radius r* where the break-off occurs and the dislocation 
spacing d versus misorientation. 



transformation 
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where rid is the number of dislocations present in the 
system. Note that this pool radius differs from the dis- 
location core radius r$ extracted from the fits to the 
Read-Shockley law. As shown in the inset of Fig. [T51 
the pool radius at the break-off point, r*, is not propor- 
tional to the dislocation spacing. Furthermore, the size 
of the liquid pools, as defined by Eq. (j4"Tj) (which is of 
course equivalent to a Gibbs construction performed for 
a cylinder around a dislocation instead of a flat homo- 
geneous liquid layer), is not uniquely determined by u, 
but also depends on the misorientation. This indicates 
that the above picture needs to be refined in order to bet- 
ter understand the condition for the coalescence of liquid 
pools above the melting point. 

As a last point, it should be recalled that diffuse- 
interface theories of grain boundaries where the grain 
orientation is treated as a scalar order parameter have 
shown the possibility that two distinct grain-boundary 
states of markedly different widths can exist at the same 
temperatur o 27 ' 28 . In contrast, aside from the dislocation- 
pairing hysteretic transition, we have found here the 
grain-boundary width to be uniquely determined at fixed 
chemical potential. However, we cannot rule out the exis- 
tence of such two-state coexistence for crystal structures 
and grain-boundary orientations other than those inves- 
tigated here, or in a narrow range of chemical potential 
very close to the melting point where numerical calcula- 
tions with the PFC model become exceedingly difficult. 
Clearly, this question warrants further investigation. 



VI. CONCLUSIONS AND OUTLOOK 

We have performed a detailed study of grain-boundary 
premelting using the phase-field crystal model. Our re- 
sults demonstrate that there is a qualitative difference be- 
tween high-angle "repulsive" and low-angle "attractive" 
grain boundaries. In the former, a continuous liquid film 
forms below the melting point, exhibiting a width that di- 
verges when the melting point is approached from below. 
For low-angle grain boundaries, melting starts at individ- 
ual dislocations. The grain boundary can be overheated 
up to a misorientation-dependent critical temperature 
at which the solid "bridges" between the liquid "pools" 
break and the system becomes liquid. Furthermore, we 
have found that a hysteretic structural transition from 
single dislocations to dislocation pairs can occur for in- 
termediate values of the misorientation. The latter, how- 
ever, is generally dependent on inclination since it is ob- 
served here for = 0° but not cj) = 30°. 

We have extracted numerically the dependence of the 
disjoining potential V(w) as a function of layer width 
w, and found that its shape is qualitatively different for 
high- and low-angle boundaries. For high angles, V(w) 
is purely repulsive for all w and reasonably well fitted 
by the exponential law of Eq. |T]), assumed in sharp- 
interface theories^ 8 -, at least for the largest misorienta- 
tion investigated here. In contrast, for low-angle grain 
boundaries, V(w) is attractive for large w, but repulsive 
for small w, and exhibits a minimum that corresponds 
to the existence of a liquid layer of finite width at the 
melting point. Furthermore, this width diverges as the 
misorientation approaches from below a critical value 9 C 
that distinguishes these two regimes. This divergence is 
smooth and reflects the progressive formation of a con- 
tinuous premelted layer by merging of liquid-like pools 
and disappearance of solid bridges between them with 
increasing misorientation. 

We have found that C is not well predicted by the ex- 
ponential form assumed in Eq. |T]) with a constant pref- 
actor. This form does not describe the large reduction of 
the grain-boundary energy due to both the decrease of 
the shear modulus at high homologous temperature, and 
local melting around dislocations that is already present 
for low- angle boundaries. We have found that, in con- 
trast, a Read-Shockley law for the grain-boundary energy 
used in conjunction with a value of the shear modulus at 
the melting point and an effective dislocation core radius, 
which describes phenomenologically dislocation-induced 
melting, yields a three to four times larger estimate of 
6 C that is in better agreement with the value obtained 
from PFC simulations. This estimate, however, is still 
too low due to the fact that dislocations are not isolated 
for 8 « 8 C as assumed in the derivation of this law. 

While this work has yielded a consistent picture of the 
thermodynamics of premelting in a microscopic model 
that can hopefully serve as a basis for developing more 
accurate mesoscopic models, it has also shown that many 
questions still need to be answered before a truly quan- 
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titative description can be obtained. First, and most 
importantly, how does the disjoining potential depend 
generally on crystal structure and grain-boundary ori- 
entation characterized by five parameters in the exten- 
sion of this work to three dimensions ? While devel- 
oping a complete theoretical description of this poten- 
tial seems difficult, there is reasonable hope that the in- 
teraction of crystal-melt interfaces due to the overlap of 
density-wave profiles for large separation (w 3> S) could 
be understood within the framework of the Ginzburg- 
Landau theor y 33 i 34 . The order parameters of this the- 
ory are the complex amplitudes of crystal density waves 
in the solid and one would expect the range of inter- 
action of the disjoining potential to be related to the 
rate of spatial decay of these density waves in the liquid. 
The fact that this theory can be derived from the PFC 
model and related quantitatively to experiments and MD 
simulations for isolated crystal-melt interfaces 3 - 4 - suggests 
that it should provide a fruitful theoretical framework 
in which to understand fundamental aspects of grain- 
boundary premelting. In particular, an asymptotic de- 
scription of the disjoining potential for large w could in 
principle shed light on the physics of the critical wetting 
angle. 

Let us finally comment on the further perspectives of 
our work. Here, we have only investigated the structural 
aspect of grain-boundary premelting. It would be inter- 
esting to study its consequences on macroscopic proper- 
ties such as the resistance to shear. In principle, shear 
can be incorporated into the PFC model by modifying 
its equations of motion^ 2 -. Since the experimental evi- 
dence for grain-boundary premelting in pure substances 
is controversial, whereas this phenomenon is well docu- 
mented in alloys, it would also be interesting to extend 
our study to this case using recently developed PFC mod- 
els for alloys^. 
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APPENDIX A: IMPLEMENTATION OF THE 
PFC MODEL 

1. Locally conserved dynamics 



where the second equality defines the linear operator L = 
(1 — e)V 2 + 2V 4 + V 6 and the nonlinear function / = 

To avoid the numerically challenging gradient terms 
in real space, the equation of motion is solved in Fourier 
space. Multiplying both sides of the equation by exp(ikx) 
and integrating over the entire volume leads to 



dti>k = L k ^ k + f k . 



(Al) 



where the Fourier modes of the density are Vfc = 
J tp exp{ikx)dx, L k = (e — l)fc 2 + 2fc 4 — k 6 is the linear 
operator in Fourier space, and f k is the Fourier transform 
of the nonlinear function /. 

Furthermore, an implicit integration scheme is used 
which allows us to use larger time steps. Instead of solv- 
ing Eq. (|Aip directly, it can be rewritten by using the 
ansatz i/} k = u(t) exp(L k t). One then obtains 

d t Tpk = Lk exp(L k t)u(t) + (d t u) exp(L k t) 
= L k exp(L k t)u(t) + f k , 

so that dtu(t) = exp(— L k t)f k . Integrating over time 
from t to t + At gives 



u(t + At) - u(t) 



t+At 



dt' eM-Lkt')f k (t') 



and with u(t) = exp(— L k t)ip k (t) in terms of tp k 
exp[-L k (t + At)]i> k (t + At) - exp(-L k t)ip k {t) 

fi +At dt'ex V (-L k t')f k (t') . 

Even if f k is not known as a function of t, it can be 
expanded in a good approximation around t' = t, leading 
to 
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The standard locally conserved dynamics for the den- 
sity field tp is given by Eq. (fl~5|) as 

dtip = (1 - e)VV + 2V 4 ^ + VV + VV 
= U + f, 



2. Non-local globally conserved dynamics 

To accelerate the search for the equilibrium solution, 
a different non-local dynamical formulation can be used 
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where the dynamics depends globally on the density field, 
as opposed to locally in the standard conserved dynamics. 
The global conservation of the order parameter has then 
to be ensured by a Lagrange multiplier. 

The conservation condition for tp is given as 



/ 



%p(x)dx — L x L y %p = 0, 



where ip = l/(L x L y ) J ip(x) is the average density. The 
free energy, including the constraint, can then be written 



T = T 



■/'■ 



tp(x)dx — L x L v ip 



where /i is the Lagrange multiplier. 
The equation of motion becomes 



d t tp = -— +/i 



8JF_ 

Sip 



= [(e-l)-2V 2 -V 4 ] tp-tp 3 +^ 
= Lip + f, 

where now L = (e - 1) - 2V 2 - V 4 and / = -tp 3 + {i. 
In Fourier space, the linear operator and the nonlinear 
function are given as It = (e — 1) + 2fc 2 — fc 4 and fk = 
—ipf, + Afci where tp 3 is the Fourier transform of tp 3 and 
jlk is that of /i. Since fj, is a constant, /}& ex (5(fc)/i. With 
this Lfc and flk, the implicit integration scheme as given 
in Eq. (|A2[) can be used. 

The Lagrange multiplier \x can be obtained from the 
condition 



= d t ^> = 
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L X Ly 
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dtip(x)dx 
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dx , 
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dx 



[{l-()ip{x)+tP 3 {x)]dx, 



since the integral over the gradients is zero for a periodic 
system. 



APPENDIX B: BOUNDARY CONDITIONS 

The two-dimensional hexagonal periodic solution given 
by Eq. (TIT))) exhibits close-packed rows of density peaks 
along the x direction and can be described by two ba- 
sis vectors a = a (1,0) and b = a (l/2, \/3/2), where 
a = 2-ir/q = 47t/a/3 is the "lattice spacing" (the spac- 
ing between density peaks). When the entire structure 
is rotated by an angle {x — > x cos <d + y sin and 



-ssinO + j/cosG), the rotated basis vectors are 



b = 



cos 8 v 
- sin 9 , 

i cos 9 



and 



v3 

2 
x/3 



sin 9 



-isin9 + ^cos9 



(Bla) 
(Bib) 



In order to exactly fit a periodic structure into the simu- 
lation box, a displacement of once the box size along the 
box axes must correspond to an integer number of steps 
along the two basis vectors, that is, we must have 



mb = 



-ia + jb 




and 



(B2a) 



(B2b) 



where (without loss of generality) < 9 < 7r/3, and n, 
to, i and j are integer numbers. The minus signs have 
been chosen by convention such that the conditions can 
be satisfied with positive integers. 

From the components of the above vector equations 
that have a zero on the right hand side, we obtain two 
conditions for the angle, 



tan 9 
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and 
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Only angles can be simulated for which four suitable in- 
tegers can be found that satisfy both conditions. Note 
that with sufficiently large integers, any angle can be ap- 
proximated to arbitrary precision. Once the four integers 
are determined, the dimensions of the simulation box are 
given by 



L y = a 



( m \ o ^ ■ o 
n cos B — to sin B 

V 2 / 2 



% ) sin + j — cos 9 

2 2 



(B4a) 
(B4b) 



For the numerical treatment, the equations have to be 
discretized. In order to accommodate both conditions for 
the system size, in general slightly different grid spacings 
have to be used along the x and y directions since the 
ratio L x /L y may be irrational and cannot be well ap- 
proximated by a single grid spacing. We always chose 
grid spacings Ax and Ay that are close to 7r/4, as in pre- 
vious PFC studie o 19 ! 20 . We have checked by repeating 
selected runs with different choices for the discretization 
that grid effects are negligible for all the results presented 
here. 

In the presence of grain boundaries, the density field 
is no longer periodic in the y direction, and Eq. (IB4b() 
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for L y does not apply. In this case, no simple condi- 
tion for L y can be given that ensures a strain-free bulk 
solid, since this would require a detailed knowledge of 
the grain-boundary structure. However, this condition 
is not as stringent as for single crystals since there is an 
additional degree of freedom: the dislocations present at 
the grain boundaries can move along the boundaries in 
response to bulk stress until a minimum of the energy 
is reached, which implies a relaxation of the bulk stress. 
Even if there is only a finite number of dislocation po- 
sitions that correspond to a local energy minimum (this 
number scales as d/a, where d is the distance between 
dislocations), as long as L y is chosen to be much larger 
than L x , the residual bulk stresses should be very weak. 
Indeed, we have varied L v by small amounts for several 
sets of parameters, and never found significant variations 
in the free-energy density. 

It should be noted that the condition on L x , Eq. (|B4aP 
still applies. For low-angle grain boundaries, the num- 
bers m and n can easily be related to explicit dislocation 
models. For instance, consider the <b = 0° inclination: 



the number m corresponds to the number of close-packed 
planes that originate at the grain boundary for each of 
the two tilted grains; the total number of edge dislocation 
is therefore equal to 2m. In turn, n indicates the num- 
ber of steps that have to be taken along a close-packed 
row before a site that is geometrically equivalent to the 
starting site can be reached by m steps along the basis 
vector b. While the average spacing d between disloca- 
tions is therefore always equal to L x /(2m), the minimum- 
energy configuration does not always correspond to equal 
spacings between dislocations. For instance, for m = 1, 
n even yields two dislocations that are evenly spaced, 
whereas n odd corresponds to a grain boundary where a 
slightly larger and smaller spacing alternate along the in- 
terface. Such configurations are well known** and consti- 
tute a local energy minimum. We did not notice any con- 
siderable difference between the behaviors of these two 
types of grain boundaries. This is to be expected since, 
due to the condition on L x in Eq. (|B4b|) . the system is 
still globally strain-free far from the grain boundary. 
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The analogous quantity in the (N, p, T) ensemble for fixed 

pressure and varying temperature is u — (T — T m )/(L/c p ), 

where L is the latent heat of fusion and c p the specific heat 

at constant pressure. 

We have G = F/[2(l + a) 



the shear modulus, bulk modulus, and Poisson's ratio in 
three dimensions, respectively. In two dimensions, we have 
<72 = c/(l — c) and Y2 = Y/(l — a 2 ). Inserting these re- 
lations in Eq. (|36[1 transforms the prefactor into 0^2/(8^), 
which is identical to the result of Ref.— , taking into ac- 
count that the modulus of the Burgers vector of a standard 
edge dislocation in the two-dimensional hexagonal struc- 
ture is equal to the lattice spacing. 



where G, Y, and a are 



